Classical Heat Transport in Anharmonic Molecular Junctions: Exact Solutions 
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We study full counting statistics for classical heat transport through anharmonic/nonlinear molecular junc- 
tions formed by interacting oscillators. Analytical result of the steady state heat flux for an overdamped an- 
harmonic junction with arbitrary temperature bias is obtained. It is found that the thermal conductance can be 
expressed in terms of temperature dependent effective force constant. The role of anharmonicity is identified. 
We also give the general formula for the second cumulant of heat in steady state, as well as the average geomet- 
ric heat flux when two temperatures are modulated adiabatically. We present an anharmonic example for which 
all cumulants for heat can be obtained exactly. For a bounded single oscillator model with mass we found that 
the cumulants are independent of the nonlinear potential. 

PACS numbers: 05.60.cd, 05.70.Ln, 44.10.+i 



Introduction — Anharmonicity or nonlinearity plays an 
important role in many physical processes. For example, it is 
crucial for umklapp scattering which gives rise to finite a ther- 
mal conductivity for bulk systems as pointed out by Peierls 
[1]. It is also an essential ingredient required for develop- 
ing functional phononic devices [2], such as thermal rectifier 
[3, 4] and heat pump [5-7]. 

Comprehensive understanding of the general behaviors of 
nonlinearity in transport is a fundamental problem in con- 
densed matter physics and nonequilibrium statistical physics. 
Although elegant approaches exist to deal with ballistic trans- 
port in harmonic systems [8, 9], solving an anharmonic model 
analytically is extremely challenging due to the intrinsic non- 
integrability of equations of motion. As a result, for classical 
heat transport, the major approaches are limited to molecular 
dynamics simulations [10, 11], perturbation theories such as 
Boltzmann equation [1, 12], mode coupling theory [13, 14], 
and/or mean-field theory [15, 16]. There does not exist an- 
alytic solutions of heat transport problem for general anhar- 
monic potentials. 

Moreover, in heat transport problems, not only the heat flux, 
but also the higher order fluctuations (noises) are of significant 
interest as they contain valuable information about the trans- 
port properties and satisfy universal symmetry known as Fluc- 
tuation Theorems (FT) [17-19]. Therefore, methods to obtain 
these higher order moments and to verify this symmetry for 
arbitrary anharmonic potentials are also highly desirable. 

In this article, we present a general approach to nonlinear 
heat transport by working in an enlarged space, which in- 
cludes the heat Q as an extra variable in the Fokker-Planck 
(FP) equation. The heat current can be expressed as the ex- 
pectation value of certain operator. If the steady state of the 
problem is known, then the analytic expression for the current 
can be obtained. We show this is the case for two different 
anharmonic oscillator systems, i.e., the overdamped interact- 
ing double-oscillator model (model I) and the bounded single- 
oscillator model (model II). The complete cumulants, as well 



as the geometric contribution to the heat flux, can be obtained 
if the full spectrum of the original FP equation is solvable. 
One such example is the V-shaped potential V(y) = k\y\. In 
addition, the Gallavotti-Cohen (GC) symmetry is verified for 
arbitrary anharmonic interactions for both models. 

We focus on the more complicated model I here and exam- 
ine the effect of nonlinearity in heat transport in this interact- 
ing system. For model II, only the results will be discussed. 

Formalism — Model I consists of two coupled Brownian 
oscillators which are connected to two Langevin heat baths 
(Fig. 1(a)). The dynamics is described by the Langevin equa- 
tions: niiXi = -d Xi V(xi - x 2 ) - 7i± l + (i = 1,2), 
where rrii and Xi are the mass and displacement of oscillator 
i. 7j is the coupling strength between the oscillator i and the 
bath i, Xi is the temperature of bath i and £j is a white noise 
with variance (£i(i)£j(t')) = 2 r y i T i 8 i j5{t — t') (k B is set to 




FIG. 1 : Schematic picture for the process of mathematical formal- 
ism for model I. (a) The original problem: two coupled oscillators 
connected to heat baths, (b) The reduced problem: single particle in 
contact with single heat bath, (c) The ultimate problem: Schrodinger 
equation for single particle. 
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1). V depicts the interaction potential, which is assumed to 
be bounded and depend only on the distance between the two 
oscillators. 

This model is a minimal model for heat transport in inter- 
acting oscillator systems. At low Reynolds number, we can 
neglect the inertia and thus obtain the overdamped dynamics 
[20]: j 1 x 1 +V'(x 1 -x 2 ) = £i and j 2 x 2 - V'(x 1 -x 2 ) = 6, 
in which the prime ' denotes the derivative. The total heat 
transferred from bath 1 to bath 2 up to time t, denoted as Q(t) 
satisfy Q = V'{x\ — x 2 )x\. 

By introducing the relative coordinate y = x\ — x 2 , we ob- 
tain two stochastic differential equations for the random vari- 
ables y and Q: 



V = 



V'(y) V'(y) 
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(fi-V'(y)). 

(1) 

Their joint probability distribution p(y, Q, t) then evolves ac- 
cording to the FP equation [21, 22] 
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where i = -i- + ± and 21 = 21 + Za. 
7 71 72 7 71 72 
By taking the Fourier transformation on Q in Eq. (2), 

we obtain the equation for the characteristic function (CF) 

z(y,ix,t) = J p(y,Q,t)e lx< ^dQ, where J is short for 

throughout the paper. Moreover, we introduce A = i\ and use 

tilde (~) to express the A-dependence (in such a notation, A is 

omitted in the parameter list). Then the equation reads 



d t z(y,t) = Lz(y,t) 



(3) 
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We shall also use A to be short for A evaluated at A 
e.g., L = L(X = 0) = L(°) and z(y,t) = J p(y,Q,t)dQh 
the probability distribution for y. Moreover, for A = 0, the 
equation dtz(y, t) = Lz(y, t) actually depicts a single Brow- 
nian oscillator connected with a single heat bath at an effective 
temperature T and subjected to a potential V (see Fig. 1(b)). 
Therefore, the stationary solution will be in the form of Boltz- 
mann distribution z(y, t) -» P^iv) = -^ e ~ V( ~ y ^ T with 
Zt = J e~ v ( v " 2 dy as the partition function. 

We are only interested in the heat Q. Therefore, by integrat- 
ing over y, we obtain the CF of Q as Z(t) = J z(y, t)dy = 
J P Q (Q,t)e x ®dQ, with P Q {Q,t) = / p(y,Q,t)dy. The cu- 
mulant generating function (CGF) reads G = In Z, which 
generates the nth-order cumulant as ((Q n )) = d^G 

A— 



Obviously, if we are able to solve Eq. (3) for general poten- 
tial V, then the problem can be completely solved, which is 
the case for harmonic potential [23]. For anharmonic cases, a 
general method is not available. 

However, using separation of variables z(y, t) = (p(y)F(t), 
Eq. (3) can still be casted into an eigenvalue problem 



F(t) 



-fit 



L(p{y) = -p<p(y). 



(5) 



with natural boundary conditions for y ±00. The struc- 
ture of the solutions is clearer if we introduce a transforma- 
tion <p(y) = e-nvVZT^y) wkh f = (21 + 2k)/(i + J_ _ 



72 ' 



-7i 



72 



—J- A). By direct substitution, it turns out that ^{y) satisfy the 
following time independent Schrodinger equation 



If d 2 
S 7 \ T dy^ + Vs 



(6) 



which describes a single particle moving in the effective po- 
tential Vs (see Fig. 1(c)), given by 



V s 



V' 2 
4T 



4^7 /1 1 A 
7i + 72 \Ti T 2 



(7) 



Then the well-known results in quantum mechanics for single 
particle can be directly evoked to discuss our situation. Two 
of them are that all the eigenvalues are real (regard A as a real 
quantity) which can be ordered as po < pi < ■ ■ ■ , and 
upon normalization, all the eigenfunctions form a complete 
orthonormal basis set, which satisfy = 5 nm and 

Er=o *n(y)*n(tf') = S(y - y'). We use (/, g) to denote the 
inner product / f{y)g(y)dy. 

The original operator L is not Hermitian, i.e., L' ^ L. 
However, it can be shown that 4> n = e y / 27 W„ satisfy U<j>„ = 
pn<t>n, namely, <j> n is an eigenfunction of L' with eigenvalue 
p n . When normalized properly, we will have 



i) = S nm , '^2<f>n(y)<Pn(y') = 8(y-y'). (8) 



n=0 



Then we can write down the formal solution of Eq. (3) as 



n=0 



(9) 



where c n is time independent and can be determined using the 
initial condition, c„ = ^0„, z(y, t = 0)^ . 

Dynamic Heat Fluxes — Based on the formalism above, 
we come to discuss the heat fluxes in nonequilibrium steady 
state. It should be pointed out that the discussions here are 
also applicable to the dynamic heat fluxes [7, 23] when system 
parameters are modulated adiabatically. 

If the long-time behavior is of main interest, then the com- 
ponent with smallest eigenvalue (po) will dominate, so that 



z(y,t) 



c e 



" Aoi ^o(y) and Z{t) -> Soe"^ ' / (po{y)dy. 
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Consequently, the CGF hxZ(t) —> —jlot and the nth cumu- 
lant of the heat fluctuations are given by 



lim 



«Q"» 



t 



(10) 



A=0 



Of course, one would not be able to get all the cumulants un- 
less fio can be obtained as a function of A. However, if one 
is only interested in the first cumulant, i.e., the average heat 
jt = — 9aAoIa=o' men nere i s an 
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flux J = lim t ^. c 
approach. 

We differentiate Eq. (5) with respect to A and obtain 

{d\L)if + Ldxfio = -(<9aMo)<A) - Mo^a^o- 



Then we perform (<fio,---) on both sides. Noticing that 

(4>o, Ld\<p ) = (i f o , d\tp ) = -jlo(<i>o, d\(po), we obtain 



J 



A=0 



4>o, (d\L)(p 



A=0 



(12) 

This equation is the analogue of the Hellmann-Feynman the- 
orem in quantum mechanics. Note that (f>o and ipo in Eq. (12) 
are evaluated at A = 0, therefore, are associated with the FP 
operator L = L^ '. It can be shown that the smallest eigen- 
value for L(L^) is /xo = and the corresponding eigenfunc- 
tions are tfaiy) oc e~ v ^l T oc Pf,* and <f>o(y) cx 1, respec- 
tively [21]. 

After normalization, Eq. (12) becomes J = 
J P^(y)dy. In this sense, serves as the heat flux 
operator, which gives the average heat flux when acted on the 
stationary state. Using L (1) = ^iT{V" - V' 2 
we finally obtain 



7i 7i 



(V" 



7i + 72 



(Ti-Ta), (13) 



where (A) T = J AP^dy is the canonical ensemble average 
of A at temperature T and we have used (V' 2 ) T = T {V") T . 
Remember that T = (72 7\ + 7iT 2 )/(7i + 72). 

Equation (13) is exact for a general potential V. It is also 
applicable to the cases with arbitrarily large temperature dif- 
ferences, i.e., the cases in a far-from-equilibrium state. As 
far as we know, this is the first exact result for anharmonic 
heat transport with general interaction, although it is only for 
a small system with overdamped dynamics. 

For harmonic potential V = \ky 2 , (V") T = k, which is 
just the force constant and is independent on the temperature. 
For anharmonic potentials, (V") T will generally be tempera- 
ture dependent. In analogy with the harmonic case, we define 
kv(T) = (V") T as the effective force constant. 

We notice that the temperature dependence of the effective 
force constant plays a key role in many of the phenomena for 
heat control and management. For example, a thermal rectifier 
is a device that directs the flow of heat [3, 4]. To use this 
oscillator system as a thermal rectifier, it is required that 



J(T! =T H ,T 2 = T L ) 



J(T 1 =T L ,T 2 = T H ) 



71 +12 

( V " ) T _ 22ZL+21I1L 



+ \. (14) 



Obviously, this condition is satisfied only when ky (T) is tem- 
perature dependent and the system is asymmetric, namely, 

71 7^ 72- 

For another example, in [24], it was demonstrated that heat 
can be shuttled across a symmetric system when the temper- 
ature of one heat bath is modulated adiabatically but with 
zero average temperature bias, e.g., T\ = T<z + f(t) with 
fit) = 0. If our system is modulated in such a way with 71 = 

72 = 27 and fit) being a periodic function of t, one can ob- 
tain an average heat flux J = j^jt- J^ p kv (^^jf^j f(t)dt 

in one period. Again, Q is nonzero only if fcy(T) de- 
pends on T. Actually, when /(t) <C T 2 , we can expand 
kvi 2T2 ^ ) into Taylor series and obtain the leading term 
as Q w k' v (T2)/ 2 (£)/87, where ' denotes derivative and 7TT 
denotes time average over one period. This pump effect is 
more significant if fcy(T2) is larger; and is vanishing when 
k' v (T 2 ) = 0. 

Since the effective force constant is independent of T for 
harmonic potential, the phenomena mentioned above cannot 
be observed using harmonic potential. On the other hand, is 
harmonic potential the only potential that has such a prop- 
erty? The answer is no. An example is the Toda potential 
V = k(e~ ax + ax) for which (V") T = ka 2 . Interestingly, 
both the harmonic lattice and Toda lattice are models display- 
ing ballistic heat transport. These facts imply a clue that tem- 
perature independent fcy (T) is related to ballistic transport. 

It is worth mentioning that there are other approximated 
mean field theories trying to convert an anharmonic potential 
in lattices to effective harmonic potential with temperature de- 
pendent force constant. For example, for the FPU-/3 model 
with V = -y-y 2 + tj-y A , the effective phonon theory (EPT) 
estimates it as fc EPT (T) = fc 2 + fc 4 ((y 4 ) T )/((y 2 ) T )[15, 25] 
while the self-consistent phonon theory (SCPT) estimates it 

as k SCPT iT) = § (k 2 + yfk 2 + 12/c 4 t)[16]. In Fig. 2, we 
compare these results with the present exact result for the 
two-oscillator model, i.e., fc oxact (T) = (fc 2 + 3fc 4 y 2 ) T . It 
is shown that both the EPT and SCPT underestimate the non- 
linear effect, although at low temperature the deviations are 
tiny. At the high temperature regime, all these theories result 
in the ~ cT 1 ' 2 dependence. For fc 2 = fc 4 = 1, the coeffi- 
cient c is 6r(3/4)/r(l/4) w 2.028 for the exact value. For 
the EPT c = r(l/4)/2r(3/4) w 1.479 and for the SCPT 
c= V3 w 1.732, which are 27% and 14% less than the exact 
value, respectively. 

Upto now we only discussed the first cumulant, i.e., the av- 
erage heat flux. However, using the standard perturbation the- 
ory, all higher order cumulants can be derived. For example, 
the second cumulant reads (see for example [26]) 



lim 



((Q 2 )) 



t 



= 2L 



(2) 



r(l) r(l) 



(15) 



where A m n = (0m, A(p n ) for short. Note that all the terms in 
the formula are again evaluated at A = 0. 
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FIG. 2: Effective force constant from different theories for the FPU- 
8 model with k-z = 1 and ki — 1. 



The excited states are involved in the expression for higher 
order cumulants. Therefore, to analytically obtain the higher 
order cumulants, we should be able to completely solve all the 
eigenstates of the Schrodinger equation Eq. (6) with A = 0, 
which describes a single particle with mass T moving in the 
potential V s = V n /4T - V"/2 (see Fig. 1(c)). This is of 
course an impossible mission for most of the anharmonic po- 
tentials. However, numerically solving it by diagonalizing the 
operator L$ is straightforward. In this sense, our method pro- 
vides a new numerical approach other than molecular dynam- 
ics to investigate the nonlinear heat transport problems. 

It can also be noticed that, even we are not able to obtain 
the CGF G = —flot as a function of A for general V, we can 
still verify the GC symmetry 



G(X) = G(A/3 - A), 



(16) 



where A/3 = \JT\ — I/T2. This result can be seen from Eqs. 
(6) and (7). Because V possess this symmetry, all its eigenval- 
ues fi n and eigenfunctions should possess this symmetry. 
However, <p n and 4>n do not follow it because the transforma- 
tion e± v / 2T contains T which is not invariant and hence GC 
symmetry is not satisfied in the transient stage. It is only valid 
in the long-time limit when the largest eigenvalue dominates. 

Geometric Heat Fluxes — It has been discovered that in 
classic heat transport, the geometric curvature might con- 
tribute to heat transport when two or more system parameters 
are modulated adiabatically [23]. Especially, when the two 
bath temperatures are modulated, there is no geometric con- 
tribution to the average heat flux for the harmonic potential. 
However, for the anharmonic FPU-/3 model, this effect was 
observed numerically. 

To speculate the role of anharmonic potential V on geomet- 
ric heat flux, we start with the expression for average heat flux 
which reads 



J geom — lim 



geom/ 



1 



Tp J Jt!T 2 ®^ 



A=0 



dT x dT 2 



(17) 

where Tt x t 2 is the curvature in the (Ti, T2) parameter space 

-^TiT 2 = (d Tl <t>o,d T2 <t>^j - (d T2 4> ,d Tl <po) . (18) 



After some derivation, we can finally obtain 

1 



T ± T 2 



A=0 



T 2 (7i +72) f-' fl Tl 
v ' m—1 



MO 



(19) 



where we still use A mn = (<fr m , A(p n ). However, since both 
V" and V are functions of y instead of operators, we can also 
evaluate it as A mn = {^ m ,A^ n ) = A nm . Now we are 
able to distinguish the effect of anharmonicity to the geomet- 
ric heat flux. For harmonic potential V = \ky 2 , V"o m = 
k (^0, fy m ), which is obviously zero for m 7^ 0. Therefore, 



A=0 



at any point (T\ , T 2 ) in the parameter space, 

and finally J ge0 m = 0. For anharmonic potentials, this term 
would generally be nonzero, which results in a non-vanishing 
curvature and geometric heat flux. 

Similar to Eq. (15), it can be seen that only when we can 
obtain all the excited states, are we able to get the final result 
of Eq. (19). However, for weak nonlinear cases, Eq. (15) and 
(19) can be used to obtain the approximated values. Consider 
the case V(y) = \k.2y 2 + all (y), with a small a. By regard- 
ing all (y) as the perturbation, we can evaluate Eq. (19) to the 
first order of a as 



T ± T 2 



A=0 



/4 0) (7i 



72)T 2 



(20) 



where (0) in the superscript indicates the unperturbed cor- 
respondences, which are for the harmonic potential. For 
the FPU-/3 model with U(y) = k 4 y 4 /4, Eq. (20) becomes 
rj-Q, which is independent of the temperature. There- 



37&4 



2( 7 l+72)fe| 

fore, the heat transfered in one period should be proportional 
to the area enclosed by the driving path. 

When one choose to modulate 7^2 while fix T\ = T 2 = T, 
it can still be proved there is no geometric heat flux. This con- 
clusion has also been found in [7] for quantum heat transport 
and in [23] for classic case with harmonic potential. 

Example — Next, we discuss an anharmonic model that 
is analytically solvable. Its potential has the V-shaped form 
V(y) = k\y\. Correspondingly, Vs = k 2 /AT — kS(y), which 
is an infinitely deep well potential plus a constant shift. This 
shift will only lift the eigenvalues and will not affect the eigen- 
states. Therefore, the full spectrum of the Schrodinger equa- 
tion with A = 0, as well as all the cumulants, should be solv- 
able. 

Actually, for this potential, we can directly solve Eq. (6). It 
is easy to show that only the ground state is bounded and the 
corresponding eigenvalue is 



Mo 



fc 2 TiT 2 
T(ji + 72) 



A 



1 



1 



X. 



(21) 



which is a quadratic form. Therefore, in the long-time limit, Q 
will be a Gaussian random variable, which only has the lowest 



two cumulants: lim t 

((Q 2 )) 
t 



^ = 7^Trm-^) = Jand 



(7l+72)T 



lim* 



2k 2 T t T 2 
T(7i+72)" 
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To obtain the geometric heat flux, the ground state eigen- 
functions are needed, which read 



~ k \y\ 



T 1 /~l 1 +T 2 /-f2 



-k\y 



2T 



l/Tl+l/72-^^l/Tl 
Tl/-i 1 +T 2 /-f 2 



(22) 

It is then easy to obtain Tt^Tq. = — 71 + 72 ^- Interestingly, 
this curvature is proportional to A. Therefore, only the first 
cumulant is nonzero. In other words, in the long-time limit, 
the geometric heat Q geom is not a random number at all. 

Model II — The idea presented here is not restricted to this 
overdamped two-oscillator model, but generally applicable to 
all heat transport problems. If the steady state can be solved, 
then at least the average heat flux can be obtained. This is 
the situation also for a bounded single-oscillator model with 
inertia, following the Langevin dynamics my = —V'(y) — 
IxV - 12V + 6+6 and Q = (-71?) + £i)y. Similar to 
Eq. (12) it can be shown that the average current is given by 



J 



7i72 



m(7i + 72) 



{Ti - T 2 ), 



(23) 



which is independent of the potential V . The immediate ques- 
tion that can be raised is whether the higher cumulants are also 
independent of VI In fact, for this model, the long-time CGF 
can be explicitly obtained as 



7i+72 
2m 




TiT 2 A 



1 



1 



(24) 

which is indeed independent of V, provided it is bounded. 
Therefore, the on-site potential has no effect on any cumulant 
of heat. Previously, this formula has been obtained for the 
harmonic case using different approaches [22, 27] and numer- 
ically tested for the FPU-/? potential [27]. 

Summary — Using the Fokker-Planck equation for the full 
counting statistics for heat, we have developed a formalism 
to obtain the dynamic and geometric heat current for arbi- 
trary anharmonic potential. The importance of anharmonic- 
ity in manipulating heat is investigated and it is found that the 
temperature dependent force constant plays a key role. The 
higher cumulants can also be obtained if the eigenspectrum of 
the original FP operator is known. As an example, we 
choose V shaped potential and obtained all the cumulants ex- 
plicitly. We also extend this method by including the inertia 
term for a bounded single particle and found that the cumu- 
lants are independent of the nonlinear onsite potential. For 
arbitrary nonlinear potential, the Gallavotti-Cohen symmetry 
is verified for both models. 

In common sense, the nonlinear transport problems are gen- 
erally thought to be not exactly solvable. We hope the exis- 



tance of exact solutions, presented here, for nonlinear trans- 
port problems will motivate further investigations on solving 
more general anharmonic systems, such as heat transport in 
nonlinear chain models which is one of the ultimate goals in 
the community of heat transport. 
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